%4-sep-2007. This code reads Julia and ACe data sets
% pick days with cont. common measurements
% 

%load c:\manoj\projects\plasma\HUA_PIU1.mat
load c:\manoj\projects\plasma\Julia_PLP.mat Julia_W fday n_w;w=n_w;
w = w - repmat(nanmean(w),[749,1]); %remove the daily mean variation
%ace_all = load('C:\Manoj\projects\ace\AceJuliaDataAndGraphs\AceJuliaDataAndGraphs\ACECombined.txt');
%load c:\manoj\projects\ace\burke_evs.mat;
load c:\manoj\projects\ace\OMNI_ELEC ace_all; % JUST ACE IEF
load c:\manoj\projects\ace\acebz ace_gsm_bz fday_gsm; % JUST ACE IEF
%load c:\manoj\projects\ace\Em.mat; %merging e field = -Vsw * Bt * sin^2(theta/2))
load c:\manoj\geomag\indices\aplist.mat;
%%

nft=64;

%Season Fliter
L = season_filter(ace_all(:,1),11,1,1,1);
 ace_all(~L,:) = [];
%season filter over---------

ace_fday = floor(ace_all(:,1));

Hd = chebi11_filter();
N_seg = 1;
N_data = 0;
mjd_date = datenum(2000,1,1);
count_juli_data  = 0;

for i = 1: 749,
  %  if sum(diff(Julia_W(i).k))+1 == length(Julia_W(i).k) &...
     if  sum(isnan(w(i,:))) <= 90, 
    L = ace_fday == Julia_W(i).fday;
    if sum(L) > 0,
     
    ace_time = ace_all(L,1)+(4*5)/(60*24);%Try to advance 20 ace time minutes
    ace_data = ace_all(L,2);

%
    L = isnan(ace_data);
    if sum(L) < 50 & sum(L) ~= 0,
        ace_data = interp1(ace_time(~L),ace_data(~L),ace_time);
    end;
%

    fprintf('Number of missing points on fday %d = %d\n', Julia_W(i).fday, sum(L));
   %if abs(mean(diff(ace_time))-0.0028) <= 2.7778e-005, % Use this with Burke data
   if abs(mean(diff(ace_time))-0.0035) <= 1e-004, %Use this with ACE min averages
        JULI = detrend(w(i,Julia_W(i).k));
        ACED = detrend(interp1(ace_time, ace_data, fday(i,Julia_W(i).k)));
        JULI1 = filter(Hd,((w(i,Julia_W(i).k))));
        ACED1 = filter(Hd,((interp1(ace_time, ace_data, fday(i,Julia_W(i).k)))));
       

        L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
        mean_ap = mean(ap(L));

        L = fday_gsm >= fday(i,Julia_W(i).k(1)) & fday_gsm <= fday(i,Julia_W(i).k(end));
        mean_imf_bz = nanmean(ace_gsm_bz(L));
         count_juli_data  = count_juli_data  + 1;
         ACED_len(count_juli_data) = length(ACED);
         JULI_len(count_juli_data) = length(JULI);
        if length(ACED) >= nft & length(JULI) >= nft & sum(isnan(ACED)) < 1,% &...
%                mean_ap >= 30,% & mean_imf_bz <= 1,
%        if sum(isnan(ACED)) < 1 & mean_ap >= 20,
%         if sum(isnan(ACED)) < 1 & mean_imf_bz <= 1 & mean_ap >=20,

            JULI_SEG(N_seg:N_seg+nft-1) = JULI(1:nft);
            ACE_SEG(N_seg:N_seg+nft-1) = ACED(1:nft);
            TIME_SEG(N_seg:N_seg+nft-1) = fday(i,Julia_W(i).k(1:nft));
            N_seg = N_seg+nft;
            N_data = N_data+1;
            data_length(N_data) = length(JULI);
        else,
           
         end;
   else,
       fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
           
    end;
    end;
end;
end;

%[Cxx,F] = mscohere(JULI_SEG,ACE_SEG,hanning(72),0,72,0.0033333);
[Cxx,F] = mscohere(JULI_SEG,ACE_SEG,hanning(nft),0,nft,0.0033333);
 
figure1 = figure;

axes('Parent',figure1,'XTick',[10 20 30 40 50 60 70 80 90 100],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);

box('on');
hold('all');
xlabel('period in minutes');
ylabel('coherence');
hold on;
semilogx((1./(60*F)),Cxx,'r');


%  PLOT SEGMENTS AND CORRELOGRAM   
% sel_data = zeros([1,N_data]);
% for i = 1:N_data,
% subplot(211);
% plotyy(TIME_SEG((i-1)*nft+1:i*nft),ACE_SEG((i-1)*nft+1:i*nft),TIME_SEG((i-1)*nft+1:i*nft),JULI_SEG((i-1)*nft+1:i*nft));
% subplot(212);
% xaxis = -(nft-1):(nft-1);
% corrl = xcorr(ACE_SEG((i-1)*nft+1:i*nft),JULI_SEG((i-1)*nft+1:i*nft));
% plot(xaxis,corrl);
% title(sprintf('%d',xaxis(find(corrl==max(corrl)))));
% grid;
% men1 = menu('','->','Exit','Sel','<-');
% if men1 == 3,
%     sel_data(i) = 1;
% end;
% end;
% 
% 
% ACE_SEG_MAT = reshape(ACE_SEG,[nft,N_data])';
% JULI_SEG_MAT = reshape(JULI_SEG,[nft,N_data])';
% TIME_SEG_MAT = reshape(TIME_SEG,[nft,N_data])';
% 
% ACE_SEG_MAT(~sel_data,:) = [];
% JULI_SEG_MAT(~sel_data,:) = [];
% TIME_SEG_MAT(~sel_data,:) = [];
% 
% ACE_SEG_NEW = reshape(ACE_SEG_MAT',[1,nft*length(ACE_SEG_MAT)]);
% JULI_SEG_NEW = reshape(JULI_SEG_MAT',[1,nft*length(JULI_SEG_MAT)]);
% TIME_SEG_NEW = reshape(TIME_SEG_MAT',[1,nft*length(TIME_SEG_MAT)]);
%     
% [Cxx1,F] = mscohere(JULI_SEG_NEW,ACE_SEG_NEW,hanning(nft),0,nft,0.0033333);
% 
% for i = 1:N_data,
% subplot(211);
% plotyy(TIME_SEG_NEW((i-1)*nft+1:i*nft),ACE_SEG_NEW((i-1)*nft+1:i*nft),TIME_SEG_NEW((i-1)*nft+1:i*nft),JULI_SEG_NEW((i-1)*nft+1:i*nft));
% subplot(212);
% xaxis = -(nft-1):(nft-1);
% corrl = xcorr(ACE_SEG_NEW((i-1)*nft+1:i*nft),JULI_SEG_NEW((i-1)*nft+1:i*nft));
% plot(xaxis,corrl);
% title(sprintf('%d',xaxis(find(corrl==max(corrl)))));
% grid;
% pause;
% end;
% 
% %save c:\manoj\projects\ace\25OCT2007
